
library (foreign)
library (survival)
library (arm)
library (pscl)
library(apsrtable)


######################################################
#
#   Negative Binomial for MID Initiation
#   Thursday, July 26 1:55 p.m.
#######################################################




## Laptop   

TermLimitsData<-  read.dta("/Users/Jeff/Dropbox/TermLimits/Replication/Data/TermLimitsAnalysisData.dta")

summary(TermLimitsData)
#TermLimitsData <- subset(TermLimitsData, year >= 1960)
#attach(TermLimitsData)


##  Model   ##

mod1<- glm.nb (countinit ~ lameduck + milservice + lameduck_mil +
                           rivalry + numbord + parliament + gender +
                                   peace + peace2 + peace3,data=TermLimitsData)
summary(mod1)
odTest(mod1)
logLik(mod1)

#Likelihood ratio test of H0: Poisson, as restricted NB model:
#Critical value of test statistic at the alpha= 0.05 level: 2.7055 
#Chi-Square Test Statistic =  4.1432 p-value = 0.0209 
#> logLik(mod1)
#'log Lik.' -800.0732 (df=12)
# So, we fail null of Poisson




forms <- strsplit(as.character(mod1$call)[2], split="~", fixed=T)
forms <- unlist(lapply(forms, function(x)paste("~", x, sep="")))
count.form <- as.formula(forms[2])

pred.dat <- data.frame(
  lameduck = c(0,1,0,1), 
  milservice = c(0,0,1,1),
  lameduck_mil = c(0,0,0,1),
  rivalry = c(0,0,0,0),
  numbord = c(2,2,2,2),
  parliament = c(0,0,0,0),
  gender = c(0,0,0,0),
  peace = c(5,5,5,5),
  peace2 = c(25,25,25,25),
  peace3 = c(125,125,125,125)
) 




X.count <- model.matrix(count.form, data=pred.dat)

#all.b <- do.call("c", mod1$coef)
set.seed(22)
b.mat <- mvrnorm(1000, mod1$coef, vcov(mod1))

count.preds <- exp(X.count %*% t(b.mat))

count.ci <- t(apply(count.preds, 1, quantile, c(.5,.025,.975)))





pred_civ_account<-mean(count.preds[1,])
pred_civ_account_low <- count.ci[1,2]
pred_civ_account_high <- count.ci[1,3]


pred_civ_tl<-mean(count.preds[2,])
pred_civ_tl_low <- count.ci[2,2]
pred_civ_tl_high <- count.ci[2,3]


pred_mil_account<-mean(count.preds[3,])
pred_mil_account_low <- count.ci[3,2]
pred_mil_account_high <- count.ci[3,3]




pred_mil_tl<-mean(count.preds[4,])
pred_mil_tl_low <- count.ci[4,2]
pred_mil_tl_high <- count.ci[4,3]


diff_civ <- count.preds[2,]-count.preds[1,]
diff_mil <- count.preds[4,]-count.preds[3,]
delta <- diff_mil-diff_civ

diff_civ_mean <- mean(diff_civ)
diff_civ_low <- quantile(diff_civ, 0.025)
diff_civ_high <- quantile(diff_civ, 0.975)

diff_mil_mean <- mean(diff_mil)
diff_mil_low <- quantile(diff_mil, 0.025)
diff_mil_high <- quantile(diff_mil, 0.975)



delta_mean<-mean(delta)
delta_low <- quantile(delta, 0.025)
delta_high<- quantile(delta, 0.975)





#########################################
## Generating Graphing Data Set
#########################################


PredictedValues <- data.frame( pred_civ_account,pred_civ_account_low,pred_civ_account_high,
                                                                         pred_civ_tl,pred_civ_tl_low,pred_civ_tl_high,
									diff_civ_mean,diff_civ_low,diff_civ_high,
                                                                         pred_mil_account,pred_mil_account_low,pred_mil_account_high,
                                                                         pred_mil_tl,pred_mil_tl_low,pred_mil_tl_high,
									diff_mil_mean,diff_mil_low,diff_mil_high,
									delta_mean,delta_low,delta_high)

write.dta(PredictedValues,file="/Users/Jeff/Dropbox/TermLimits/Replication/Models/Military/Conditional/MultipleDisputes/PredictedValues.dta")



